library(EBImage)

Take one photo as an example to explain the process of data preprocessing.

Step 0: Plot original photo with points

indices = 920
setwd("~/Documents/GitHub/Spring2020-Project3-group8/xjx")
train_dir <- "../data/train_set/"
train_image_dir <- paste(train_dir, "images/", sep="")
train_pt_dir <- paste(train_dir,  "points/", sep="")
train_label_path <- paste(train_dir, "label.csv", sep="") 
image.path_sub <- paste0(train_image_dir, sprintf("%04d", indices), ".jpg")
image.list_sub = list(EBImage::readImage(image.path_sub))
img = Image(image.list_sub[[1]], colormode = 'Color')
display(img, method = 'raster')

display(img, method="raster")
load("../output/fiducial_pt_list.RData")
fiducial_pt_list_sub <- fiducial_pt_list[indices]
add_point <- function(n){text(x = fiducial_pt_list_sub[[1]][n,1],
                              y = fiducial_pt_list_sub[[1]][n,2],
                              label = as.character(n), col = "white", cex = 0.8)}
#lapply(1:78,add_point)
for(i in 1:78){
  add_point(i)
}

Step 1: Determine symmetry axis and rotate photo

points = data.frame(x = fiducial_pt_list_sub[[1]][1:78,1], y = fiducial_pt_list_sub[[1]][1:78,2])
single = c(35,36,37,38,44,52,56,59,62)
double = data.frame(x1 = c(1,2,3,4,5,6,7,8,9,19,20,21,22,23,24,25,26,39,40,41,42,43,50,51,57,58,63),
                    x2 = c(10,15,14,13,12,11,18,17,16,31,30,29,28,27,34,33,32,49,48,47,46,45,54,53,55,60,61))
pos = data.frame(x = points[single,1], y = points[single,2])
for(i in 1:nrow(double)){
  pos = rbind(pos, data.frame(x = (points[double[i,1],1]+points[double[i,2],1])/2, 
                              y = (points[double[i,1],2]+points[double[i,2],2])/2))
}
display(img, method="raster")
points(pos$x, pos$y, col = 'red')

# Can't use y~x because coeff of x is too big.
model = lm(x~y, data = pos)
summary(model)
## 
## Call:
## lm(formula = x ~ y, data = pos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0524 -2.8936 -0.2338  2.6352 11.8967 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 612.536324   3.501275  174.95   <2e-16 ***
## y            -0.128857   0.008427  -15.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.259 on 34 degrees of freedom
## Multiple R-squared:  0.873,  Adjusted R-squared:  0.8693 
## F-statistic: 233.8 on 1 and 34 DF,  p-value: < 2.2e-16
y0 = 0:100*10
xhat = model$coefficients[1] + model$coefficients[2]*y0
display(img, method="raster")
points(xhat, y0, col = 'red')
points(pos$x, pos$y, col = 'blue')

a = 750
b = 1000
angle = atan(model$coefficients[2])
img_rotate = rotate(img, angle/2/pi*360)
display(img_rotate, method="raster")
if(angle > 0){
  points(c(a*sin(angle), 0, a*sin(angle)+b*cos(angle), b*cos(angle)),
       c(0, a*cos(angle), b*sin(angle), a*cos(angle)+b*sin(angle)), 
       col = c('green','blue','white','orange'), cex = 10)
}else{
  points(c(0, a*sin(-angle), b*cos(angle), b*cos(angle)+a*sin(-angle)),
       c(b*sin(-angle), a*cos(-angle)+b*sin(-angle), 0, a*cos(angle)), 
       col = c('green','blue','white','orange'), cex = 10)
}

if(angle > 0){
  new_points = data.frame(x = NULL, y = NULL)
  for(i in 1:78){
    x = points[i,1] - 0
    y = points[i,2] - 0
    alpha = atan(y/x)
    d = sqrt(x^2+y^2)
    newx = d*cos(alpha + angle) + 0 + 750*sin(angle)
    newy = d*sin(alpha + angle) + 0 
    
    new_points = rbind(new_points, data.frame(x = newx, y = newy))
  }
}else{
  new_points = data.frame(x = NULL, y = NULL)
  for(i in 1:78){
    x = points[i,1] - 0
    y = points[i,2] - 0
    alpha = atan(y/x)
    d = sqrt(x^2+y^2)
    newx = d*cos(alpha + angle) + 0
    newy = d*sin(alpha + angle) + 0 + 1000 * sin(-angle)
    
    new_points = rbind(new_points, data.frame(x = newx, y = newy))
  }
}
display(img_rotate, method="raster")
new_add_point <- function(n){text(x = new_points[n,1],
                              y = new_points[n,2],
                              label = as.character(n), col = "white", cex = 0.8)}
#lapply(1:78,new_add_point)
for(i in 1:78){
  new_add_point(i)
}

new_pos = data.frame(x = new_points[single,1], y = new_points[single,2])
for(i in 1:nrow(double)){
  new_pos = rbind(new_pos, data.frame(
    x = (new_points[double[i,1],1]+new_points[double[i,2],1])/2, 
    y = (new_points[double[i,1],2]+new_points[double[i,2],2])/2))
}
display(img_rotate, method="raster")
points(new_pos$x, new_pos$y, col = 'red')

# Can't use y~x because coeff of x is too big.
new_model = lm(x~y, data = new_pos)
summary(new_model)
## 
## Call:
## lm(formula = x ~ y, data = new_pos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0142 -2.8701 -0.2087  2.6191 11.7634 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 607.651852   3.875266 156.803   <2e-16 ***
## y            -0.000301   0.008289  -0.036    0.971    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.224 on 34 degrees of freedom
## Multiple R-squared:  3.879e-05,  Adjusted R-squared:  -0.02937 
## F-statistic: 0.001319 on 1 and 34 DF,  p-value: 0.9712
new_y0 = 0:100*10
new_xhat = new_model$coefficients[1] + new_model$coefficients[2]*new_y0
display(img_rotate, method="raster")
points(new_xhat, new_y0, col = 'red')
points(new_pos$x, new_pos$y, col = 'blue')

Step 2: Determine the zoom ratio and scale the photo

fixdist = -170
my_fixdist = apply(new_points[c(1,10),],2,mean)[2] - new_points[44,][2]
rate = fixdist/my_fixdist
new2_points = new_points*unlist(rate)
img_zoom = resize(img_rotate, nrow(img_rotate)*unlist(rate))
display(img_zoom, method="raster")
new_add_point2 <- function(n){text(x = new2_points[n,1],
                              y = new2_points[n,2],
                              label = as.character(n), col = "white", cex = 0.8)}
#lapply(1:78,new_add_point2)
for(i in 1:78){
  new_add_point2(i)
}

Step 3: Determine translation distance and translation

img_translate = translate(img_zoom, c(500-new2_points[37,1], 375-new2_points[37,2]))
display(img_translate, method="raster")
points(500, 375, col = 'red', cex = 1, pch = 16)

new3_points = data.frame(x = NULL, y = NULL)
for(i in 1:78){
  new3x = new2_points[i,1] + (500-new2_points[37,1])
  new3y = new2_points[i,2] + (375-new2_points[37,2])
  new3_points = rbind(new3_points, data.frame(x = new3x, y = new3y))
}
display(img_translate, method="raster")
new3_add_point <- function(n){text(x = new3_points[n,1],
                              y = new3_points[n,2],
                              label = as.character(n), col = "white", cex = 0.8)}
#lapply(1:78,new3_add_point)
for(i in 1:78){
  new3_add_point(i)
}

img_final = img_translate[60:940, 60:720, ]
display(img_final, method="raster")
points(new3_points$x-60, new3_points$y-60, col = 'blue', cex = 0.8, pch = 16)
points(500-60, 375-60, col = 'red', cex = 2, pch = 16)